7.6.1. Métodos locales

Los métodos locales hacen referencia a hacer cálculos entre celdas que se sobreponen en una pila de rásters. Es decir la localización espacial de las celdas es la misma. La figura muestra cómo podríamos calcular la media de dos rásters que se sobreponen.

Hay varias funciones para hacer cálculos de este tipo. Por ahora utilizaremos la función app() para demostrar este método.

Veamos el código en R para calcular la media de la pila de tres capas que hemos creado; es decir la funcion app() en terra toma los tres valores de cada pixel de cada capa que se sobreponen y calcula la media (mean). Esto se hace en toda la extensión del área que definimos.

r_stack_media <- app(x = r_stack, # Este es el objeto que representa la pila de rásters
                     fun = mean   # Esta es la función para calcular la media 
                     )
r_stack_media

## class       : SpatRaster 
## dimensions  : 10, 10, 1  (nrow, ncol, nlyr)
## resolution  : 36, 18  (x, y)
## extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 
## source(s)   : memory
## name        : mean 
## min value   :  1.1 
## max value   : 73.7

Note que el nuevo objeto r_stack_media representa la media de las tres capas en r_stack. R asigna el nombre de la función a la capa final en esta caso name: es mean.

Visualicemos el promedio de esto tres rásters.

plot(r_stack_media)
text(r_stack_media, digits=1, cex=0.8)

Usted ya puede ver las potenciales aplicaciones de estas funciones en un pila de rásters. Por ejemplo, calcular la media o la varianza de la temperatura durante un año (una pila de 12 capas) en su área de estudio.

Ahora veamos un ejemplo utilizando datos globales reales!!!!

tmean <- geodata::worldclim_global(
    var="tavg", # temperatura media
    res=10,     # resolución gruesa
    path = "C:/Talleres_R/IntroR_espacial" # folder donde desea guardar los datos
                      )

tmean

## class       : SpatRaster 
## dimensions  : 1080, 2160, 12  (nrow, ncol, nlyr)
## resolution  : 0.1666667, 0.1666667  (x, y)
## extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326) 
## sources     : wc2.1_10m_tavg_01.tif  
##               wc2.1_10m_tavg_02.tif  
##               wc2.1_10m_tavg_03.tif  
##               ... and 9 more source(s)
## names       : wc2.1~vg_01, wc2.1~vg_02, wc2.1~vg_03, wc2.1~vg_04, wc2.1~vg_05, wc2.1~vg_06, ... 
## min values  :    -45.8840,   -44.80000,   -57.92575,   -64.19250,   -64.81150,   -64.35825, ... 
## max values  :     34.0095,    32.82425,    32.90950,    34.19375,    36.25325,    38.35550, ...

par(mfrow=c(1,2)) 

terra::plot(tmean$wc2.1_10m_tavg_01, main="Temperatura media - Enero", range=c(-50, 40), mar=c(4, 3, 4, 3), plg=list(shrink=0.9, cex=.8), 
    pax=list(side=1:2, cex.axis=.6))
terra::plot(tmean$wc2.1_10m_tavg_10, main="Temperatura media - Agosto", range=c(-50, 40), mar=c(4, 2, 4, 4), legend=FALSE, 
    pax=list(side=c(1,4), cex.axis=.6))

tmean_sd <- stdev(x = tmean, # SpatRaster object
                  pop = T, # Si es verdadero (T), entonces la desviación estándard de la población es calculada
                  na.rm = T # Si es verdadero (T), entonces valores NA son ignorados
                  )

plot(tmean_sd)

tmean_mean <- terra::app(x = tmean,  # SpatRaster
                  fun = mean) # función
plot(tmean_mean, main="Temperatura media anual (°C)")

7.6.2. Métodos zonales y globales

Los métodos zonales y globales hacen referencia a hacer cálculos del área total del ráster utilizando la función global() o zonas de un ráster utilizando la función zonal() (estas zonas se pueden definir con otro capa). Por ejemplo, la figura muestra cómo podemos calcular la sumatoria de una zona en el ráster usando un polígono.

Veamos un ejemplo de cómo podemos utilizar la base de datos global de la huella humana (Human footprint, en inglés) para calcular el promedio de huella humana para cada país del mundo.

La huella humana reune la mayor parte de las acciones humanas con el potencial de deteriorar la naturaleza (ver el trabajo de Venter et al 2016 para mayores detalles).

La escala es 0 = muy baja, hasta 50 = muy alta huella humana.

Note que en nuestros códigos vamos a utilizar el simbolo %>%, conocido como forward-pipe operator en inglés, que envia el resultado de una expresión a la otra función dentro de un mismo objeto. Este operador%>% fue originalmente desarrollado paquete magrittr.

En este caso cargamos la base de datos “footprint” y después la agregamos 10 veces utilizando la función aggregate(). Tenga en cuenta que la función aggregate() se encuentra en varios paquetes, entonces es buena práctica especificar que estámos utilizando la función del paquete terra, asi: terra:aggregate().

¿Por qué agregamos el ráster huella humana? Porque queremos la huella humana en un ráster de resolución más gruesa para acelerar el procesamiento de nuestros siguientes pasos (es decir menos tiempo de procesamiento). Más adelante veremos unos ejemplos del tiempo de procesamiento en diferentes resoluciones.

huella <- geodata::footprint(year = 2009,
                             path ="C:/Talleres_R/IntroR_espacial")%>% 
          terra::aggregate(fact=10, fun=mean)

## 
|---------|---------|---------|---------|
=========================================
                                          

huella

## class       : SpatRaster 
## dimensions  : 2160, 4320, 1  (nrow, ncol, nlyr)
## resolution  : 0.08333333, 0.08333333  (x, y)
## extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326) 
## source(s)   : memory
## name        : wildareas-v3-2009-human-footprint 
## min value   :                              0.00 
## max value   :                             45.21

Ahora carguemos los países del mundo. Es una base de datos que contiene información vectorial. Más adelante dedicaremos una sección a este formato.

paises <- geodata::world(resolution = 5,
                         level = 0,
                        path = "C:/Talleres_R/IntroR_espacial")
paises

##  class       : SpatVector 
##  geometry    : polygons 
##  dimensions  : 231, 2  (geometries, attributes)
##  extent      : -180, 180, -90, 83.65625  (xmin, xmax, ymin, ymax)
##  coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
##  names       : GID_0      NAME_0
##  type        : <chr>       <chr>
##  values      :   ABW       Aruba
##                  AFG Afghanistan
##                  AGO      Angola

Ahora si, apliquemos la función zonal() del paquete terra para calcular el promedio de la huella humana en cada país del mundo.

huella_paises <- terra::zonal(huella,
                              paises,
                              as.raster=T, # si es F obtiene solo los datos.
                              fun = mean
                              )
plot(huella_paises)

7.6.3. Seleccionado un área geográfica dentro de un ráster

Es muy común que usted tenga un ráster con una extensión global, pero usted quiera enfocar su análisis en una región específica. Para esto vamos a utilizar la función crop() del paquete terra::.

Por ejemplo, seleccionemos los datos de la huella humana solo para Colombia.

Por ahora transformemos la información vectorial a ráster, utilizando la función rasterize() del paquete terra::.

Note que vamos a utilizar el ráster huella como plantilla para rasterize() los polígonos paises; es decir obtendremos un ráster de las mismas dimensions:, resolution:, extent y coord. ref.:.

Para los valores del ráster utilizaremos el nombre de la columna en el objeto vectorial que contiene el nombre de los países, es decir NAME_0.

paises_r <- terra::rasterize(paises, huella, field="NAME_0")
paises_r

## class       : SpatRaster 
## dimensions  : 2160, 4320, 1  (nrow, ncol, nlyr)
## resolution  : 0.08333333, 0.08333333  (x, y)
## extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326) 
## source(s)   : memory
## categories  : NAME_0 
## name        :      NAME_0 
## min value   : Afghanistan 
## max value   :       Åland

plot(paises_r)

Las categorias del objeto paises_r pueden ser inspeccionadas con la función levels() del paquete terra. Generalmente son dos columnas, la primera es un numero entero que identifica los valores de la celda y la segunda es una etiqueta que indica el nombre de la categoria. Si revisa la tabla podrá ver que Colombia es el valor 43.

knitr::kable(head(terra::levels(paises_r)))%>%
  kableExtra::kable_paper("hover", full_width = F) %>%
  kableExtra::scroll_box(width = "100%", height = "200px")

Lo que vamos hacer es copiar el objeto paises_r y llamarlo Colombia, luego vamos a reclasificar todos los valores diferentes a 43 con NA (sin valores), es decir escoger Colombia.

Colombia <- paises_r
Colombia[Colombia != 43] <- NA
Colombia

## class       : SpatRaster 
## dimensions  : 2160, 4320, 1  (nrow, ncol, nlyr)
## resolution  : 0.08333333, 0.08333333  (x, y)
## extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326) 
## source(s)   : memory
## categories  : NAME_0 
## name        :   NAME_0 
## min value   : Colombia 
## max value   : Colombia

# Compare el código de arriba con este de abajo. ¿Cuales son las diferencias? ¿Podría ser de utilidad?

# Colombia <- paises_r == 43

Visualicemos el mapa de Colombia. Hemos seleccionado lo que queremos, pero el extent: sigue siendo el mismo, es decir todo el globo entero, por eso vemos el país tan pequeño.

plot(Colombia)

Podemos utilizar la función trim() para remover las filas y columnas con valores NA (alrededor de los valores del mapa de Colombia), y como consecuencia se cambia el extent: del nuevo objeto.

Colombia_t <- terra::trim(Colombia)
Colombia_t

## class       : SpatRaster 
## dimensions  : 200, 146, 1  (nrow, ncol, nlyr)
## resolution  : 0.08333333, 0.08333333  (x, y)
## extent      : -79, -66.83333, -4.25, 12.41667  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326) 
## source(s)   : memory
## categories  : NAME_0 
## name        :   NAME_0 
## min value   : Colombia 
## max value   : Colombia

plot(Colombia_t)

Ahora sí utilicemos la función crop() para extraer los datos de temperatura media en el extent: del objeto Colombia_t.

Colombia_huella <- terra::crop(huella, Colombia_t)
Colombia_huella

## class       : SpatRaster 
## dimensions  : 200, 146, 1  (nrow, ncol, nlyr)
## resolution  : 0.08333333, 0.08333333  (x, y)
## extent      : -79, -66.83333, -4.25, 12.41667  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326) 
## source(s)   : memory
## name        : wildareas-v3-2009-human-footprint 
## min value   :                              0.00 
## max value   :                             37.86

Visualicemos el nuevo objeto. Como podemos ver la función crop() extrae lo valores de la huella humana con base en el extent: del objeto Colombia_t.

plot(Colombia_huella)

Para tener solamente los valores de la huella humana dentro de los limites de Colombia_t debemos utilizar la función mask(). Veamos:

Colombia_huella_mask <- terra::mask(Colombia_huella, Colombia_t)
plot(Colombia_huella_mask)

Ahora bien, por qué no utilizar la opción mask() desde un principio? Hay varias razones.

• las dos rásters deben ser del mismo extent:. Al intentar implementar la función mask() directamente al mapa global huella obtenemos el Error. Es decir el mapa huella (es global) y el mapa Colombia_t tienen diferente extent:.

Colombia_huella_2 <- terra::mask(huella, Colombia_t)

## Error: [mask] extents do not match

• la segunda razón es que es más eficiente para R implementar la función crop(), porque selecciona filas y columnas del área definida por el extent: que queremos. Al contrario, mask() require hacer cálculos más complejos para seleccionar los pixeles que están bordeando la figura del mapa de Colombia.

7.6.4. Seleccionando valores en un ráster

En ocasiones necesitamos utilizar solo unos valores determinados dentro del ráster.

Por ejemplo, queremos obtener un mapa con las temperaturas por encima de 25 grados. Para lo cual podemos reclasificar los valores menores a 25 con NA (no values).

tmean_caliente <- tmean_mean
tmean_caliente[tmean_caliente < 25] <- NA
tmean_caliente

## class       : SpatRaster 
## dimensions  : 1080, 2160, 1  (nrow, ncol, nlyr)
## resolution  : 0.1666667, 0.1666667  (x, y)
## extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326) 
## source(s)   : memory
## name        :    mean 
## min value   : 25.0000 
## max value   : 30.9879

plot(tmean_caliente)

Pero tambien podemos necesitar un rango , por ejemplo entre 15 y 20 grados. Para eso utilizamos el operador | que indica “o” (los valores menores a 15 “o” valores mayores a 20, serán iguales a NA).

tmean_templado <- tmean_mean
tmean_templado[tmean_templado < 15 | tmean_templado >20] <- NA
plot(tmean_templado)

Tambien podemos extraer lo valores en una pila de rásters. Por ejemplo, seleccionemos los valore mayores a 10 grados en todas las rásters.

tmean_frio <- tmean
tmean_frio[tmean_frio > 10] <- NA
tmean_frio

## class       : SpatRaster 
## dimensions  : 1080, 2160, 12  (nrow, ncol, nlyr)
## resolution  : 0.1666667, 0.1666667  (x, y)
## extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326) 
## source(s)   : memory
## names       : wc2.1~vg_01, wc2.1~vg_02, wc2.1~vg_03, wc2.1~vg_04, wc2.1~vg_05, wc2.1~vg_06, ... 
## min values  :     -45.884,       -44.8,   -57.92575,    -64.1925,    -64.8115,   -64.35825, ... 
## max values  :      10.000,        10.0,    10.00000,     10.0000,     10.0000,    10.00000, ...

7.6.5. Seleccionando valores en un ráster con otro ráster

Vamos a extraer los valores de precipitación global en areas con temperaturas mayores a 20°C (tmean_caliente).

Primero, descarguemos la precipitación utilizando el paquete geogata:

prec <- geodata::worldclim_global(
          var="prec",
          res=10,
          path = "C:/Talleres_R/IntroR_espacial"
                                 )
prec

## class       : SpatRaster 
## dimensions  : 1080, 2160, 12  (nrow, ncol, nlyr)
## resolution  : 0.1666667, 0.1666667  (x, y)
## extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326) 
## sources     : wc2.1_10m_prec_01.tif  
##               wc2.1_10m_prec_02.tif  
##               wc2.1_10m_prec_03.tif  
##               ... and 9 more source(s)
## names       : wc2.1~ec_01, wc2.1~ec_02, wc2.1~ec_03, wc2.1~ec_04, wc2.1~ec_05, wc2.1~ec_06, ... 
## min values  :           0,           0,           0,           0,           0,           0, ... 
## max values  :         908,         793,         720,        1004,        2068,        2210, ...

Ahora calculemos la precipitación media anual:

prec_mean <-terra::app(prec, mean)
prec_mean

## class       : SpatRaster 
## dimensions  : 1080, 2160, 1  (nrow, ncol, nlyr)
## resolution  : 0.1666667, 0.1666667  (x, y)
## extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326) 
## source(s)   : memory
## name        :     mean 
## min value   :   0.0000 
## max value   : 932.5833

terra::plot(prec_mean)

Finalmente, utilicemos la función mask() para extraer los valores de precipitación en areas con temperaturas mayores a 20°C.

prec_caliente <- terra::mask(prec_mean, tmean_caliente)
prec_caliente

## class       : SpatRaster 
## dimensions  : 1080, 2160, 1  (nrow, ncol, nlyr)
## resolution  : 0.1666667, 0.1666667  (x, y)
## extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326) 
## source(s)   : memory
## name        :     mean 
## min value   :   0.0000 
## max value   : 932.5833

Visualizemos los dos rásters:

# Creemos dos columnas para los plots
par(mfrow=c(1,2))

# Temperature > 20 °C
plot(tmean_mean*0, legend=F) # Fondo gris para mostrar los continentes 
t_breaks <- seq(minmax(tmean_caliente)[1], minmax(tmean_caliente)[2], 1) # escala de la temperatura
terra::plot(tmean_caliente, add=T, type="interval", plg=list(x="bottom", cex=1.2, title="°C"), breaks = t_breaks)
title("Temperatura media > 20 °C", line = -2, cex.main=1.6)

# Precipitación en áreas con temp  > 20 °C
p_breaks <- seq(minmax(prec_caliente)[1], minmax(prec_caliente)[2], 150)
plot(tmean_mean*0, legend=F)
terra::plot(prec_caliente,  add=T, type="interval", plg=list(x="bottom", cex=1.2, title="mm"), col=rev(topo.colors(15)), breaks=p_breaks)
title("Precipitación en áreas con temp  > 20 °C", line = -2, cex.main=1.6)

7.6.6. Extraer los valores de un ráster

En algún momento de sus análisis necesitará los valores de los rásters para hacer algunos cálculos adicionales o preparar gráficas. El paquete terra ofrece la función values() para extraer lo valores del ráster.

Veamos un ejemplo con la temperatura media y la precipitación global.

Visualicemos una gráfica que muestre los valores pareados de temperatura media y precipitación.

Interesante ver que las precipitaciones medias anuales en áreas frías son menores que áreas más calientes.

plot(values(tmean_mean), values(prec_mean))

Ahora tomemos solamente los datos de precipitación para Colombia.

Colombia_t_prec <- terra::crop(prec_mean,  # precipitación media
                        terra::aggregate(Colombia_t, 2) # agreguemos colombia 2 veces para que sea la misma resolución de la precipitación.
                        )%>%
                  resample(aggregate(Colombia_t, 2))%>% # asegurarnos que las rásters sean iguales en extent, resolution, etc.
                  mask(aggregate(Colombia_t, 2))       # Seleccionar solamente el borde administrativo de Colombia

plot(Colombia_t_prec)

El Pacífico siempre lluvioso!

Seleccionemos también la temperatura media para Colombia.

Colombia_t_tmean <- terra::crop(tmean_mean, aggregate(Colombia_t, 2))%>%
                    resample(aggregate(Colombia_t, 2))%>%
                    mask(aggregate(Colombia_t, 2))
plot(Colombia_t_tmean)

Ahora visualicemos los valores pareados de temperatura media y precipitación para Colombia con relación al global.

Definitivamente tropical!

  plot(values(tmean_mean), values(prec_mean), xlim=c(-55, 30), ylim=c(0, 1000),
       ylab="",yaxt="n", xlab="",xaxt="n", pch=21, col="black")

par(new=T)
plot(values(Colombia_t_tmean), values(Colombia_t_prec), col="red", xlim=c(-55, 30), ylim=c(0, 1000), pch=22)

legend("topleft", legend=c("Colombia", "Global"), col=c("red", "black"), pch=c(22,21), cex=0.8)

7.6.7. Proyectando rásters

Hasta ahora hemos utilizado coordenadas geográficas para estos análisis. Sin embargo, los análisis más locales o regionales requieren de ráster proyectadas. Es decir los datos proyectados en una superficie plana. En especial si queremos hacer cálculos de áreas para informar sobre aspectos del área de estudio (ej., área de bosques fragmentados, etc).

Hay diferentes proyecciones cartográficas, pero la que más vamos a utilizar en el análisis espacial es la que mantienen el sentido del área. Cada país y región tiene proyecciones apropiadas que preservan el area, entonces hay que hacer el ejercicio de buscar la mejor proyección para el área de estudio.

Por fortuna hay muchos sitios para buscar y spatialrefence y epsg.io.

Por ejemplo, si buscamos “Colombia” en spatialreference nos aparecerá algo asi:

En nuestros códigos podemos utilizar cualquiera de los dos formatos ESPG (es el que aparece en el primer pantallazo de ‘spatialreference’) o Proj4.

Por ejemplo si utilizamos EPSG:3116 el equivalente Proj4 es: +proj=tmerc +lat_0=4.5962004166 +lon_0=-74.077507916 +k=1 +x_0=1000000 +y_0=1000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs

Veamos cómo utilizamos esto en R.

• Carguemos la cobertura vegetal (solamente la clase árboles, “trees”),
• seleccionemos los datos para Colombia y
• utilicemos la proyección EPSG:3116 para reproyectar el ráster.

La función project() de terra incluye varios argumentos:

• y = hace referencia a “coordinate reference system, (por sus siglas en inglés), en este caso EPSG:3116.
• method = “bilinear” es util para variables continuas y “near” para variables categóricas. Vamos a utilizar “bilinear porque es la proporción de bosque es cada celda (es decir que es una variable continua).
• res = hace referencia al tamaño de celda (en metros) que queremos al final. Es decir 1Km por 1Km. Celdas cuadradas, lo que facilitará y mejorará la exactitud (accuracy, en inglés) de los cálculos.

arboles <- geodata::landcover(var= 'trees',
                    path = "C:/Talleres_R/IntroR_espacial"                            )%>%
  terra::crop(disagg(Colombia_t,2))%>%
  terra::resample(terra::disagg(Colombia_t,2))%>%
  terra::mask(disagg(Colombia_t,2))%>%
  terra::project(y="EPSG:3116", method="bilinear", res = 1000)
  
arboles

## class       : SpatRaster 
## dimensions  : 1857, 1358, 1  (nrow, ncol, nlyr)
## resolution  : 1000, 1000  (x, y)
## extent      : 451351.3, 1809351, 19203.33, 1876203  (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS / Colombia Bogota zone (EPSG:3116) 
## source(s)   : memory
## name        : trees 
## min value   :     0 
## max value   :     1

plot(arboles)

Utilicemos la función cellSize() de terra para calcular el área de cada celda de este ráster. Casi todas las celdas son cercanas a 1Km2 (~1,000m x ~1,000m), pero hay algunas muy pequeñas variaciones, en especial en el oriente de Colombia, porque esta proyección esta centrada más en la mitad de Colombia.

Continuemos con nuestro análisis, asumiendo que las celdas son de igual tamaño. Sin embargo, tenga en cuenta reportar la exactitud en sus resultados.

arboles_cellsize <- cellSize(arboles, unit="km")
plot(arboles_cellsize)

7.6.8. Reclasificando rásters

La reclasificación de un ráster es de mucha utilidad en el análisis espacial. Por ejemplo si usted quiere reclasificar la cobertura de árboles en Colombia en cuatro categorias.

Para esto debemos crear una matriz que contenga el valor inicial, el valor final del ráster y que valor queremos al final. Repetir el proceso hasta terminar el rango.

Por ejemplo, definamos tres nuevas categorías de acuerdo a la proporción de árboles y reclasifiquemolas con 1, 2 y 3. En la primera linea se dice, de 0 a 0.25 reclasifiquelo a 1.

m <- c(0, 0.25, 1,
       0.25, 0.5, 2,
       0.5, 1, 3)
rclmat <- matrix(m, ncol=3, byrow=TRUE)
rclmat

##      [,1] [,2] [,3]
## [1,] 0.00 0.25    1
## [2,] 0.25 0.50    2
## [3,] 0.50 1.00    3

Ahora apliquemos esta matriz al ráster.

rc <- classify(arboles,
               rclmat)
plot(rc, legend = FALSE,
     col = c("dark green",  "red", "yellow"), axes = FALSE)
legend("topright",
       legend = c("Alta", "Media",  "Baja"),
       fill = c("dark green", "yellow",  "red"),
       border = FALSE,
       bty = "n") # turn off legend border

Podríamos tener unas datos del área de cada categoría de la siguiente forma:

terra::freq(rc)

##   layer value  count
## 1     1     0     91
## 2     1     1 213491
## 3     1     2 164580
## 4     1     3 761372

Resolución espacial

Creemos 3 grillas de la misma extensión pero diferente tamaño de celda:

pixel <- c(1, 10, 100, 1000, 5000, 7500, 10000) # Tamaño de celda

pixel_size <- list()
for(i in 1:length(pixel)){
  
pixel_size[[i]] <- rast(ncol=pixel[i], nrow=pixel[i])
  values(pixel_size[[i]]) <- 1:ncell(pixel_size[[i]])

}

Ahora apliquemos una función sencilla a cada ráster y calculemos el tiempo que tarda cada una en procesar.

tiempo_proceso <- list()
tiempo <- list()

for(j in 1:length(pixel_size)){
  

  tiempo_proceso[[j]] <- system.time(sqrt(pixel_size[[j]]))[3]
    
}

## 
|---------|---------|---------|---------|
=========================================
                                          

|---------|---------|---------|---------|
=========================================
                                          

tiempo <- unlist(tiempo_proceso)
 
df1 <- data.frame(pixel, tiempo)

df1

##   pixel tiempo
## 1     1   0.00
## 2    10   0.00
## 3   100   0.00
## 4  1000   0.01
## 5  5000   0.22
## 6  7500   2.17
## 7 10000   3.92

Visualicemos en una gráfica para que pueda observar la tendencia.

plot(df1$pixel,df1$tiempo, xlab = "Tamaño de celda", ylab = "Tiempo de procesamiento", )

Por ahora tengamos en cuenta que el tiempo de procesamiento es una limitante cuando manejamos rásters. Las soluciones a este problema van desde la optimización de sus códigos hasta la utlización de supercomputadores para acelerar los procesos.